function [price,outRate,inRate]=snowball_MC2(S0,C,knockinB,knockoutB,r,Sigma,T,Ntraj,Nstep)
%%%%%%%%%%%%%函数说明%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%传统蒙特卡罗方法计算雪球产品价格%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%观测频率说明：雪球敲出每月观测；敲入每日观测；
%%%%%%%%%%输入参数：
%%%%%%%%%%S0：标的初始价格；knockinB：敲入界限；knockoutB:敲出界限；C:年化票息率；
%%%%%%%%%%r：无风险利率；Sigma：标的波动率；T：到期时间；
%%%%%%%%%%Ntraj：模拟路径数；Nstep：每条路径一年的时间分割数
%%%%%%%%%%输出参数：
%%%%%%%%%%Price:雪球产品价格;outRate:敲出路径所占比例;inRate:敲入未敲出路径所占比例
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%EXAMPLE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<1
T=1;S0=1; r=0.03; Sigma=0.13; knockinB=0.85;knockoutB=1.03;C=0.2;
Ntraj=1000; Nstep=364;
end
%%%%%%%%%%%%%%%%%%%%%%%筛选knock out的观测日%%%%%%%%%%%%%%%%%%%%
% 标记每个月最后一天，并且生成一个365*1的vector, 只有在每个月的最后一天
% Vector取1，否则取0；
startDate = datetime('2023-01-01', 'Locale', 'zh_CN');
endDate = datetime('2023-12-31', 'Locale', 'zh_CN');
dateVector = datevec(startDate:endDate);
monthDays = accumarray(dateVector(:,2), dateVector(:,3), [], @max);
idx = zeros(12,1);
for i = 1:12
    idx(i) = find(dateVector(:,2)==i&dateVector(:,3)==monthDays(i));
end
Vector = zeros(365, 1);
Vector(idx) = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S=simuPath(S0,T,r,Sigma,Ntraj,Nstep);%一列为一条价格路径
payoff=zeros(1,Ntraj);
holdingPeriod=T*ones(1,Ntraj);%记录雪球产品的有效持有期

%Case I: knock out
outORnot=S.*Vector>S0*knockoutB;%检验节点价格是否出现敲出
outPath=find(any(outORnot,1));
period=[0,1:T*Nstep]'.*outORnot(:,outPath);%用于寻找敲出路径的有效持有期
period(period==0)=NaN;
holdingPeriod(outPath)=min(period,[],'omitnan')/Nstep;
notoutPath=setdiff([1:Ntraj],outPath);%未敲出
% Case II: knock in but not out
inORnot=S<=S0*knockinB; % filter knock in
inPath=find(any(inORnot,1));
innotoutPath=setdiff(inPath,intersect(inPath,outPath));% filter in but not out
% case III: filter not in not out
notinnotoutPath=setdiff([1:Ntraj],[inPath,outPath]);%未敲入未敲出

%%%%%%%计算每一种情况的收益
payoff(outPath)=S0*(C*holdingPeriod(outPath));
payoff(notinnotoutPath)=S0*(C*holdingPeriod(notinnotoutPath));
payoff(innotoutPath)= S0*min(S(Nstep,innotoutPath)./S0-1,0);

price=exp(-r*T)*mean(payoff);%计算所有路径收益均值并进行贴现
outRate=length(outPath)/Ntraj;%敲出的比率
inRate=length(inPath)/Ntraj;%敲入的比率

% visualization
disp(['nums of out path is ', num2str(length(outPath))]);
disp(['nums of knock in but not out path is ', num2str(length(innotoutPath))])
disp(['nums of not in not out path is ', num2str(length(notinnotoutPath))])

% 最小的敲出时间
%[row,col,v] = find(outORnot);
%first_one_in_row = accumarray(col, row, [], @min);
% 使用find函数找到矩阵中值为true的元素的位置信息
% 平均敲出时间
%knock_out_time_Avg = sum(first_one_in_row)/(30*length(outPath));
end